#pragma once


void make_ft_spincolor(char *momentum_order,spincolor *S,spincolor *FT, int *xsource,  int sign)
{

  if (strcmp (momentum_order,"TXYZ") == 0)
  {
        printf("TXYZ momentum order \n");
  for(int imom=0;imom<Nmom;imom++) memset(FT[imom],0,sizeof(spincolor));

  double lP[4];
  int x[4];
 for(int imom=0;imom<Nmom;imom++)
 {
   for(int idir=0;idir<4;idir++) lP[idir]=P[imom][idir];

   for(x[0]=0;x[0]<T;x[0]++)
    for(x[1]=0;x[1]<L;x[1]++)
      for(x[2]=0;x[2]<L;x[2]++)
        for(x[3]=0;x[3]<L;x[3]++)
        {
          int dest=x[3]+L*(x[2]+L*(x[1]+L*x[0]));
          double arg=0;
          for(int idir=0;idir<4;idir++) arg+=lP[idir]*(x[idir]-xsource[idir]);
          complex eipx=cos(arg)+I*sign*sin(arg);
          for (int is=0; is<4 ; is++)
                for (int ic=0; ic<3; ic++)
			 FT[imom][is][ic]=FT[imom][is][ic]+S[dest][is][ic]*eipx;

             }
       }

    }
// En principio como FT=sum_x S e^ipx deberia poder usar la misma definición en los dos casos
// pues P viene como una variable externa que en cada caso se ha definido con el orden de momentos que toca
// Pero mantengo la diferencia para posibles futuras pruebas

  if (strcmp (momentum_order,"TYXZ") == 0)
  {
        printf("TYXZ momentum order \n");
          for(int imom=0;imom<Nmom;imom++) memset(FT[imom],0,sizeof(spincolor));

  double lP[4];
  int x[4];
 for(int imom=0;imom<Nmom;imom++)
 {
   for(int idir=0;idir<4;idir++) lP[idir]=P[imom][idir];

   for(x[0]=0;x[0]<T;x[0]++)
      for(x[2]=0;x[2]<L;x[2]++)
        for(x[1]=0;x[1]<L;x[1]++)
          for(x[3]=0;x[3]<L;x[3]++)
        {
          int dest=x[3]+L*(x[2]+L*(x[1]+L*x[0]));
          double arg=0;
          for(int idir=0;idir<4;idir++) arg+=lP[idir]*(x[idir]-xsource[idir]);
          complex eipx=cos(arg)+I*sign*sin(arg);
          for (int is=0; is<4 ; is++)
                for (int ic=0; ic<3; ic++)
                         FT[imom][is][ic]=FT[imom][is][ic]+S[dest][is][ic]*eipx;

             }
       }
  }

  if (strcmp (momentum_order,"TYZX") == 0)
  {
        printf("TYZX momentum order \n");
	  for(int imom=0;imom<Nmom;imom++) memset(FT[imom],0,sizeof(spincolor));

  double lP[4];
  int x[4];
 for(int imom=0;imom<Nmom;imom++)
 {
   for(int idir=0;idir<4;idir++) lP[idir]=P[imom][idir];

   for(x[0]=0;x[0]<T;x[0]++)
      for(x[2]=0;x[2]<L;x[2]++)
        for(x[3]=0;x[3]<L;x[3]++)
          for(x[1]=0;x[1]<L;x[1]++)
        {
          int dest=x[3]+L*(x[2]+L*(x[1]+L*x[0]));
          double arg=0;
          for(int idir=0;idir<4;idir++) arg+=lP[idir]*(x[idir]-xsource[idir]);
          complex eipx=cos(arg)+I*sign*sin(arg);
          for (int is=0; is<4 ; is++)
                for (int ic=0; ic<3; ic++)
                         FT[imom][is][ic]=FT[imom][is][ic]+S[dest][is][ic]*eipx;

             }
       }
  }
}



void make_ft_colorcolorspinspin(char *momentum_order, colorcolorspinspin *S,colorcolorspinspin *SFT, int *xsource,  int sign)
{
  for (int is1=0; is1<4; is1++)  for (int ic1=0; ic1<3; ic1++)
  {
    int sc=3*is1+ic1;
    // Temporal spincolor
    spincolor *s=(spincolor*)malloc(sizeof(spincolor)*(T*L*L*L));

    printf("====> FT is=%d, ic=%d, isc=%d\n",is1, ic1, sc);

    for (int t=0; t<T; t++)    for (int x=0; x<L; x++)   for (int y=0; y<L; y++)     for (int z=0; z<L; z++)
    {
          int dest=z+L*(y+L*(x+L*t));
          for (int is2=0; is2<4; is2++) for (int ic2=0; ic2<3; ic2++)
                 s[dest][is2][ic2]=S[dest][ic2][ic1][is2][is1];
     }

  // temporal spincolor to save the fourier transform of each external spin-color
  spincolor *sFT_tmp=(spincolor*)malloc(sizeof(spincolor)*(Nmom));
  // call the main fourier transform routine
  make_ft_spincolor(momentum_order,s,sFT_tmp,xsource, sign);
  // copy to colorcolorspinspin
  for (int imom=0; imom<Nmom; imom++)
        for (int is2=0; is2<4; is2++)for (int ic2=0; ic2<3; ic2++)
                SFT[imom][ic2][ic1][is2][is1]=sFT_tmp[imom][is2][ic2];

  free(s); free(sFT_tmp);
 }
}
